%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                                 Model                                   %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ---------------------------- Description -------------------------------- 

    % This script computes the main objects of the model given a vector of
    % parameters. It does essentially the same thing as the
    % Calibration_Function. However, outside of the function environment,
    % it allows us to store the objects of the model more easily.

    % In addition, it computes some objects of the model that are not used
    % for calibration but are important for the analysis, such as the
    % entrant employment distribution and the exit rates by firm age.

%--------------------------------------------------------------------------
%                            Preliminaries                                %
%--------------------------------------------------------------------------

    fprintf('------------------------ Started Computation of Model ------------------------\n')
    fprintf('\n\n')

    % Efficiency parameter for an entrant's product creation.
    p.phie = paramguess(1);

    % Step in productivity of a variety replaced by creative destruction  
    p.l = paramguess(2);                             

    % Efficiency parameter for an incumbent's product creation.
    p.phix = paramguess(3);

    % Moment of the exogenous productivity distribution raised to (sigma-1)
    p.wbar = paramguess(4);

    % Exogenous product destruction rate
    p.d0 = paramguess(5);

    % ---------------------------------------------------------------------
    %                 1. Solving for Direct Parameters                    %
    % ---------------------------------------------------------------------

    % For an overview of this part of the code, please refer to Section 2.2
    % of the Replication Package

    %----------------------------------------------------------------------

    % Rate of incumbent product creation - see Proposition (M-3)
    p.xstar = (p.phix/(p.zeta*p.phie))^(1/(p.zeta-1));

    % Solves for alpha and I  - see Section (RP-2.2)
    Direct_Parameters;

    % Creative destruction rate - see Proposition (M-3)
    p.tau = p.alpha*(p.eta+p.d0)/(1-p.alpha);

    % Flow entry rate - see Proposition (M-3)
    p.z = 1/p.alpha*p.tau-p.xstar;

    % Net rate of product accumulation - see Section (RP-2.2.2)
    p.psi = p.xstar - (p.tau+p.d0);

    % Average Efficiency gain of new products raised to sigma-1 - see Equation (M-5)
    p.qbar = p.alpha*p.l^(p.s-1)+(1-p.alpha)*p.wbar;

    % The CES markup level - see Section (M-2)
    p.mu = p.s/(p.s-1);

    %------------------------- Computing Moments --------------------------

    % Compute the Aggregate Productivity Growth rate - see Equation (RP-1)
    g = p.eta/(p.s-1) + p.In + p.tau/(p.alpha*(p.s-1))*(p.qbar-1);

    % Computes the numerator of Firm Age Density - see Equation (RP-4)
    Firms_function = @(a) p.psi*exp((p.psi-p.eta).*a)./(p.psi - p.xstar*(1-exp(p.psi.*a)));

    % Integrates the numerator over firm age for denominator - see Equation (RP-4)
    Firms = integral(Firms_function, 0, Inf);

    % Computes firm age density as a function - see Equation (RP-4)
    FirmAge_dens = @(a) Firms_function(a) / Firms;

    % Computes the average number of products per firm - see Equation (RP-5)
    xProds = (p.z * Firms)^(-1);

    % Computes the aggregate entry rate - see Equation (RP-6)
    EntryRate = p.z * xProds;

    %---------------------- Product Distribution --------------------------

    % Preallocate a vector for the unconditional product distribution
    ProductDist = zeros(1, p.maxprods);

    % Compute the Gamma function - see Equation (RP-3)
    GammaFunc = @(a) p.xstar*(1-exp(p.psi*a))./(p.xstar*(1-exp(p.psi*a))-p.psi);

    % Computes the unconditional product probabilities recursively...
    for n = 1:p.maxprods

        % ... by defining the integrand - see Equations (RP-2) and (RP-7)
        Integrand = @(a) (1-GammaFunc(a)) .* GammaFunc(a).^(n-1) .* FirmAge_dens(a);

        % ... and summing over firm age - see Equations (RP-7)
        ProductDist(n) = integral(Integrand, 0, Inf);

    end

    % Clear the integrand function
    clear Integrand

    % Compute the Exit Rate implied by the model - see Equation (RP-52)
    ExitRate = (p.tau + p.d0) * ProductDist(1);

    % Share of firms older than 10 years - see Equation (RP-56)
    Firms_10y = integral(FirmAge_dens, 10, Inf);

    % ---------------------------------------------------------------------
    %                  2. Solving for Average Firm Size                   %
    % ---------------------------------------------------------------------  

    % For an overview of this part of the code, please refer to Section 2.3
    % of the Replication Package

    %----------------------------------------------------------------------

    % The constant C shows up in the differential equation for k(Delta) - see Section (RP-2.3.1)
    ConstC = p.rho + p.d0 + (p.eta+p.d0)*(p.qbar/(1-p.alpha)-1);

    % Boundary condition delivers k(s/s-1) - see Section (RP-2.3.1)
    k_func_mu = 1/(ConstC*(p.s-1))*(p.s/(p.s-1))^(-p.s);

    % Constant of integration for the diff. equation - see Equation (RP-12)
    ODE_const = (1/(ConstC*(p.s-1))-p.s/(ConstC*(p.s-1)+p.In*(p.s-1)^2)+1/(ConstC+p.In*p.s))*(p.s/(p.s-1))^(-p.s-ConstC/p.In);

    % Solution of differential equation to obtain k(lambda) - see Equation (RP-11)
    k_func_lam = p.l^(-p.s)*(p.l/(ConstC+p.In*(p.s-1))-1/(ConstC+p.In*p.s)) + ODE_const*p.l^(ConstC/p.In);
    
    % -------------- Joint Distribution of Gaps and Productivity ----------

    % Scale parameter for Pareto distribution - see Equation (RP-16)
    p.wmin = (p.wbar*(1-(p.s-1)/p.varrho))^(1/(p.s-1)); 

    % Growth rate of average efficiency - see Section (RP-2.3.2)
    p.gQ = (p.tau/p.alpha)*(p.qbar-1)/(p.s-1) + p.In;

    % Computes joint distributions and misallocation terms - see Section (RP-2.3.2) 
    [FNC, F1, M, Lam, Deltagrid, p.qgrid_adj] = Joint_Dist(p);

    % Compute ls^P (times Misallocation terms) - see Equation (RP-9)
    LPNwMLam = (1 - (p.zeta-1)*p.xstar/(p.zeta*(p.rho+p.tau+p.d0)))/(p.phie*(p.alpha*k_func_lam*p.l^(p.s-1)+(1-p.alpha)*k_func_mu*p.wbar));

    % Compute ls^P by removing the misallocation terms accordingly
    LPN = LPNwMLam * M^(p.s-1)*Lam^(p.s);

    % With ls^P, we can compute l = L/N - see Equation (RP-8)
    LFN = LPN + ((p.eta+p.d0)/(1-p.alpha) - (p.zeta-1)*p.xstar/p.zeta)/p.phie;

    % Having ls^P and l, compute production share s^P on the BGP
    prod_share = LPN/LFN;

    % Compute the Average Firm Size - see Equation (RP-22)
    AvFirmSize = LFN * xProds;

    % ---------------------------------------------------------------------
    %                 3. Mark-up and Sales growth                         %
    % ---------------------------------------------------------------------

    % For an overview of this part of the code, please refer to Section 2.4
    % of the Replication Package

    %----------------------------------------------------------------------

    % Conditional distribution of product age given firm age and products - see Section (RP-2.4.4)
    [AgeDist, avec, Pmatrix, p0] = Conditional_Dist(p);

    % Take differences to obtain the associated density
    TimePeriods = p.maxTime/p.periodlength + 1;
    AgeDensity = zeros(TimePeriods, TimePeriods, p.maxprods);
    AgeDensity(2:end,:,:) = diff(AgeDist, 1, 1);
    AgeDensity(1,:,:) = AgeDist(1,:,:);

    % Conditions conditional product probabilities on firm survival
    CondProbs =  Pmatrix ./ repmat((1-p0), p.maxprods, 1);

    % ------------------------ a. Mark-ups --------------------------------

    % Markup by product age for varieties with competitors - see Equation (RP-24)
    markupbyage_C = min(p.l * exp(p.In*avec),(p.s/(p.s-1)));

    % Expected markup by product age - see Equation (RP-25)
    markupbyage = p.alpha*markupbyage_C + (1-p.alpha)*p.mu;

    % Expected markup given firm age and number of products - see Equation (RP-23)
    Emarkup_An = permute(sum(AgeDensity .* repmat(markupbyage', 1, TimePeriods, p.maxprods), 1),[3 2 1]);
    
    % Expected markup given firm age - see Equation (RP-23)
    Emarkup = sum(Emarkup_An .* CondProbs, 1);

    % Finds the index on age grid associated with a firm of 10 years
    [~, ind10] = min(abs(avec-10));

    % Computes markup growth by age 10 - see Equation (RP-30)
    MarkupGrowth = Emarkup(ind10) - Emarkup(1);

    % ------------------------ b. Sales -----------------------------------

    % Expected product-level sales by product age - see Equation (RP-26)
    salesbyage = exp((p.s-1)*(p.In-p.gQ).*avec) / (M^(p.s-1)*Lam^(p.s-1)) .* (p.alpha*p.l^(p.s-1)*markupbyage_C.^(1-p.s) + (1-p.alpha)*p.wbar*p.mu^(1-p.s));

    % Expected product-level sales given firm age and number of products - see Equation (RP-23)
    Sales_An = permute(sum(AgeDensity .* repmat(salesbyage', 1, TimePeriods, p.maxprods), 1),[3 2 1]);

    % Vector with number of products to turn product level into firm level
    prodmat = repmat((1:p.maxprods)',1,TimePeriods);

    % Expected firm-level sales given firm age - see Equation (RP-23)
    Esales = sum(Sales_An .* CondProbs .* prodmat, 1);

    % Computes sales growth as a function of firm age
    Esales = log(Esales) - log(Esales(1));

    % Computes firm-level sales growth by age 10 - see Equation (RP-30)
    SalesGrowth = Esales(ind10);

    % ------------------ c. Cost Weighted Markup --------------------------

    % Expected cost-weighted markup by product age - see Equation (RP-27)
    costmarkupsbyage = (p.alpha*p.l^(p.s-1)*markupbyage_C.^(1-p.s) + (1-p.alpha)*p.wbar*p.mu^(1-p.s)) ./ (p.alpha*p.l^(p.s-1)*markupbyage_C.^(-p.s) + (1-p.alpha)*p.wbar*p.mu^(-p.s));

    % Expected cost-weighted markup given firm age and number of products - see Equation (RP-23)
    CostMarkup_An = permute(sum(AgeDensity .* repmat(costmarkupsbyage', 1, TimePeriods, p.maxprods), 1),[3 2 1]);

    % Expected cost-weighted markup given firm age - see Equation (RP-23)
    CostMarkup_A = sum(CostMarkup_An .* CondProbs, 1);

    % Integrate over firm age to get aggregate markup - see Equation (RP-30)
    Integrand = @(a) interp1(avec, CostMarkup_A, a, "linear", CostMarkup_A(end)) .* FirmAge_dens(a);
    CostMarkup = integral(Integrand, 0, Inf);

    % Clear the remaining variables
    clear markupbyage_C markupbyage Emarkup_An ind10 salesbyage Sales_An prodmat costbyage markupcostbyage costmarkupsbyage CostMarkup_An Integrand

    % ---------------------------------------------------------------------
    %                   4. Employment Distribution                        %
    % ---------------------------------------------------------------------

    % For an overview of this part of the code, please refer to Section 2.5
    % of the Replication Package

    %----------------------------------------------------------------------

    % Create a matrix with Firm Age density with expanded dimensions
    FirmAge_matrix = repmat(FirmAge_dens(avec), TimePeriods, 1, p.maxprods);

    % Obtain conditional distribution of product age given number of products - see Section (RP-2.5)
    AgeDist_N = permute(trapz(avec, AgeDist .* FirmAge_matrix, 2), [1, 3, 2]);

    % Take differences from distribution to get associated density
    AgeDensity_N = [AgeDist_N(1, :); diff(AgeDist_N, 1, 1)];

    % Computes the function D_n^1(y) - see Equation (RP-34)
    Dy1 = Employ_Dist(p, F1, AgeDensity_N, Deltagrid, LPNwMLam);

    % Exogenous employment grid over which we want employment
    Dy_grid = exp(p.yvec)';

    % Creates a matrix for D_n^n(y) - see Equation (RP-31)
    Dy_n = zeros(p.ypoints, p.maxprods);

    % For each product we consider...
    for n = 1:p.maxprods

        % ... transform distribution according to Equation (RP-31)
        Dy_n(:,n) = interp1(n*Dy_grid, Dy1(:,n), Dy_grid);

    end

    % Remove NaN associated with off-grid interpolation
    Dy_n(isnan(Dy_n)) = 0;
    
    % Aggregate the distribution over the number of products - see Equation (RP-32)
    Employment = sum(Dy_n .* repmat(ProductDist, p.ypoints, 1), 2);

    % Make sure the resulting distribution adds up to 1 by correcting numerical imprecisions
    Employment = Employment/Employment(end);
    
    % Obtain the associated density for log(employment)
    Emp_density = [0; diff(Employment)] / (p.yvec(2)-p.yvec(1));

    % Transform density of log(employment) to density of employment
    Emp_density = Emp_density ./ Dy_grid;

    % ---------------------- Employment Shares ----------------------------

    % Determine the endpoints for the binned distribution
    emppoints = [4 9 19 49 99 249 499 999 2499 4999 9999];

    % Preallocate vector to store the share of employment in each bin
    empshare = zeros(1, length(emppoints)+1);

    % Denote previous employment point
    empold = 0;

    % For each endpoint we target...
    for i = 1:length(emppoints)
 
        % Integrand is simply the density-weighted employment
        Integrand = @(y) y.*interp1(Dy_grid, Emp_density, y, 'linear', 0);

        % Compute share of employment attributable to current bin
        empshare(i) = integral(Integrand, empold, emppoints(i)) / integral(Integrand, 0, Inf);

        % Update starting point for the next bin to be analyzed
        empold = emppoints(i);

    end

    % Computes the employment share of large firms - see Equation (RP-33)
    empshare(end) = integral(Integrand, empold, Inf) / integral(Integrand, 0, Inf);
    clear empold

    % Computes the employment share of large firms - see Equation (RP-33)
    ShareLarge = empshare(end);

    % ---------------------------------------------------------------------
    %                 5. Collecting the Results                           %
    % ---------------------------------------------------------------------

    % Create a table To display calibrated parameters and model moments
    parameter_names = categorical({'phi_e'; 'lambda'; 'I'; 'delta'; 'phi_x'; 'alpha'; 'omega_bar'});
    parameter_values = round([p.phie; p.l; p.In; p.d0; p.phix; p.alpha; p.wbar^(1/(p.s-1))], 3);
    moments_names = categorical({'Avg. Firm Size';'Productivity Growth (%)';'Markup Growth (%)';'Entry Rate (%)'; 'Sales Growth (%)'; 'Weighted Average Markup (%)'; 'Emp. Share of Large Firms'});
    moments_model = round([AvFirmSize; g*100; MarkupGrowth*100; EntryRate*100; SalesGrowth*100; (CostMarkup-1)*100; ShareLarge*100], 2);
    moments_target = round([p.AverageFirmSize_Moment; p.g_Moment*100; p.MarkupGrowth_Moment*100; p.EntryRate_Moment*100; p.SalesGrowth_Moment*100; (p.Markup_Moment-1)*100; p.ShareLarge_Moment*100] ,2);
    columnnames = {'Parameter','Calibration','Moment', 'Model', 'Target'};
    calibration_table = table(parameter_names, parameter_values, moments_names, moments_model, moments_target, 'VariableNames', columnnames);
    disp(calibration_table)

    % Clear the variables used to construct the table
    vars = ({'parameter_names', 'parameter_values', 'moments_names', 'moments_model', 'moments_target', 'columnnames', 'calibration_table'});
    clear(vars{:})
    clear vars;

    % ---------------------------------------------------------------------
    %              6. Entrant Employment Distribution                     %
    % ---------------------------------------------------------------------

    % For an overview of this part of the code, please refer to Section
    % 2.5.2 of the Replication Package

    %----------------------------------------------------------------------

    % Obtain the entrant employment distribution - see Section (RP-2.5.2)
    EntrantDist = Entrant_Dist(p, F1, Deltagrid, LPNwMLam);

    % Interpolate the Entrant Distribution to find 10th quantile
    Ent_Func = @(x) interp1(p.yvec, EntrantDist, x) - 0.1;

    % Solve for the 10th quantile
    Emp_10perc = fzero(Ent_Func, 0);

    % Interpolate the Entrant Distribution to find 90th quantile
    Ent_Func = @(x) interp1(p.yvec, EntrantDist, x) - 0.9;

    % Solve for the 90th quantile
    Emp_90perc = fzero(Ent_Func, 0);
    clear Ent_Func

    % Compute the implied 90-10 ratio for entrant employment
    Ratio90_10 = exp(Emp_90perc) / exp(Emp_10perc);

    % ---------------------------------------------------------------------
    %                       7. Exit Rates by age                           %
    % ---------------------------------------------------------------------   

    % For an overview of this part of the code, please refer to Section 2.7
    % of the Replication Package

    %----------------------------------------------------------------------

    % Preallocate vector for the probability of exit given number of products
    exitrates_prod = zeros(1, p.maxprods);

    % Grid for the number of products that could be destroyed in a period
    Kprods = 1:p.maxprods;

    % For every possible number of products that a firm may have...
    for n = 1:p.maxprods

        % Vector that stores the terms to be summed over m
        Finalterm = zeros(1, p.maxprods + 1);

        % For every possible number of products could be created
        for m = 0:p.maxprods

            % Compute the term of the outer sum - see Equation (RP-40)
            outer = exp(-p.xstar*n)*(p.xstar*n)^m/factorial(m);

            % Compute the term of inner sum - see Equation (RP-40)
            inner = exp(-(p.tau+p.d0)*n)*((p.tau+p.d0)*n).^Kprods./factorial(Kprods);

            % Add the terms of the inner sum - sum starts at n+m - see Equation (RP-40)
            inner = sum(inner(n+m:end));

            % Multiply the outer term with the inner sum - see Equation (RP-40)
            Finalterm(m+1) = outer*inner;

        end

        % Stores probability of exit for a firm with n products - see Equation (RP-40)
        exitrates_prod(n) = sum(Finalterm);

    end

    % Transforms exit conditional on number of products into exit conditional on firm age - see Equation (RP-41)
    Exit_age = repmat(exitrates_prod', 1, TimePeriods);
    Exit_age = sum(Pmatrix .* Exit_age, 1);
    clear inner outer Finalterm Kprods m;

    % ---------------------------------------------------------------------
    %                                                                     %
    % ---------------------------------------------------------------------    

    fprintf('------------------------- Ended Computation of Model ------------------------\n')
    fprintf('\n\n')




    





